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ABSTRACT 

We study static neutron stars with poloidal magnetic fields and a simple class of 
electric current distributions consistent with the requirement of stationarity. For this 
class of electric current distributions, we find that magnetic fields are too large for static 
configurations to exist when the magnetic force pushes a sufficient amount of mass off- 
center that the gravitational force points outward near the origin in the equatorial plane. 
(In our coordinates an outward gravitational force corresponds to din gu/dr > 0, where 
t and r are respectively time and radial coordinates and gu is coefficient of dt'^ in 
the line element.) For the equations of state (EOSs) employed in previous work, we 
obtain configurations of higher mass than had been reported; we also present results 
with more recent EOSs. For all EOSs studied, we find that the maximum mass among 
these static configurations with magnetic fields is noticeably larger than the maximum 
mass attainable by uniform rotation, and that for fixed values of baryon number the 
maximum mass configurations are all characterized by an off-center density maximum. 

Subject headings: stars: neutron — stars: magnetic fields 



1. INTRODUCTION 

Over the years, the typical magnitudes of the surface magnetic fields of pulsars — as inferred 
from measured spindown rates and simple magnetic dipole models — have been ~ 10^^ — 10^^ G 
(Taylor, Manchester & Lyne 1993). Assuming magnetic fiux conservation, fields of ~ 10^^ G would 
naturally arise from typical main sequence star surface magnetic fields of ~ 100 G during a decrease 
in radius by a factor of ~ 10^ (Shapiro & Teukolsky 1983). At the extreme end of fields attainable 
by fiux conservation, the largerst observed white dwarf magnetic field of 5 x 10^ G leads to a neutron 
star field of ~ 1 x 10^^ G (Carroll &: Ostlie 1996), while the largest observed main sequence stellar 
magnetic field of 3.4 x 10^ G (Bohm-Vitense 1989) also suggests a possible field of a few times 
10^4 G. 
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Several independent circumstantial arguments link the class of objects known as "soft 7-ray 
repeaters" (SGRs), and perhaps the so-called "anomalous X-ray pulsars" (AXPs), with neutron 
stars having magnetic fields > 10"*^^ G — the so-called "magnetars" (Duncan & Thompson 1992; Usov 
1992; Paczyiiski 1992; Thompson & Duncan 1995, 1996; Vasisht & Gotthelf 1997). (Table 1 displays 
some observed properties of these objects.) In addition to the circumstantial arguments, more direct 
evidence for magnetic fields of 2 — 8 x 10^^ G is available for two of the five known SGRs, in the form 
of measured periods and spindown rates of associated X-ray pulsars (Kouvcliotou ct al. 1998, 1999)."^ 
Furthermore, the observed X-ray luminosities of the AXPs may require a field strength B > 10^^ G 
(Chatterjee, Hcrnquist, & Narayan 2000; Heyl & Kulkarni 1998). The population statistics of SGRs 
suggest that magnetars may constitute a significant fraction (> 10%) of the neutron star population 
(Kouveliotou et al. 1994, 1998). As mentioned above, there are isolated examples of progenitor stars 
which could yield fields of ~ 10^^ G by flux conservation, but these isolated examples do not seem 
sufficient to account for a significant fraction of the neutron star population. Thus, it seems likely 
that some mechanism generates magnetic fields in nascent neutron stars. For example, Duncan 
& Thompson (1992) suggest that the smoothing out of differential rotation and convection could 
generate fields as large as 3 x 10^''(Pi/l ms)^^ G, where Pi is the initial rotation period of the 
neutron star. 

These considerations motivate study of the effects of ultra-strong magnetic fields on neutron 
star properties. In this we have been inspired by the pioneering work of Bocquet, Bonazzola, 
Gourgoulhon &; Novak (1995), who performed relativistic calculations of axisymmetric neutron 
star structure in which the standard stress-energy tensors of a perfect fluid and the electromagnetic 
field were employed, and were comparable in magnitude. The maximum fields they found were of 
order 10^^ G, with increases of 13 to 29% in the maximum mass of nonrotating stars for various 
equations of state. 

An additional motivation is provided by the recent findings that magnetic fields of strengths 
larger than 10^^ G affect the EOS of dense matter directly through drastic changes in the com- 
position of matter (Chakrabarty 1996; Chakrabarty, Bandyopadhyay, &; Pal 1997; Yuan Sz Zhang 
1999; Broderick, Prakash, & Lattimer 2000). The EOS is altered by both the Landau quantization 
of the charged particles (such as protons, electrons, etc.) and the interactions of the magnetic mo- 
ments, including the anomalous magentic moments of the neutral particles (such as the neutron, 
strangeness-bearing A-hyperon etc.) with the magnetic field. In this work, we consider only the 
effects of the magnetic field on the structure, through its influence on the metric, in order to facil- 
itate a comparison with the earlier work of Bocquet et al. (1995). The additional changes caused 
by the direct effects of the magnetic field on the EOS will be reported in a future work (Cardall, 
Broderick, Prakash, &; Lattimer 2000). 



^For ongoing discussion of this interpretation of X-ray timing data see e.g. Marsden et al. (1999a); Woods et 
al. (1999a); Harding et al. (1999); Marsden, Rothschild, & Lingenfelter (1999b); Chatterjee, Hernquist, & Narayan 
(2000). 
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While Bocquet et al. (1995) also presented solutions for rotating neutron stars endowed with 

large magnetic fields, in this work we present only static solutions. In terms of the potential 
observability of the effects of large magnetic fields, the most relevant situation appears to be for non- 
rotating stars. Neutron stars with the highest inferred magnetic fields — the so-called "magnetars" — 
are all observed to be rotating very slowly, with periods on the order of seconds. The effects of 
such slow rotation should have a negligible impact on the neutron star structure. 

In this paper we extend the work of Bocquet et al. (1995). The theoretical formalism is outlined 
in §2, which serves to put the problem in context. In §3 we shed light on an issue left somewhat 
unclear by Bocquet et al. (1995): What physically determines the maximum mass and magnetic 
field for a given maximum density or given baryon mass? In order to explore these questions we 
have chosen to solve the structural equations using a Green's function technique rather than the 
spectral technique employed by Bocquet et al. (1995), and we also searched for the maximum mass, 
for a given magnetic field distribution, in a different way. An appendix describes our numerical 
methods and tests of our code. In §4 we present an illuminating view of constant baryon mass and 
constant magnetic moment sequences, and present higher mass configurations than those found by 
Bocquet et al. (1995) for the equations of state (EOSs) they employed. In addition, we present 
the results of analogous calculations using three more recent EOSs. Summary and outlook are 
contained in S5. 



We consider stationary neutron star models in which the equation of state is independent of 
the magnetic field. The relevant equations and some properties of neutron stars in this limit were 
studied by Bonazzola et al. (1993) and Bocquet et al. (1995). The stress-energy tensor is given by 
the sum of the standard stress-energy tensors of a perfect fluid and the electromagnetic field: 



where e and p are respectively the rest-frame energy density and pressure, u°' is the fluid 4- velocity, 
gai3 are the metric components, and F^p = ^/3,a ~ ^a,/3) where is the electromagnetic potential 
1-form. 

At least two relativistic formulations of the problem of strongly nonspherical axisymmetric stars 
have been developed. Most authors studying rapid rotation have adopted the approach of Bardeen 
& Wagoner (1971), which explicitly assumes an isotropic stress tensor and is thus incompatible with 
electromagnetic fields. Building on earlier work of Bonazzola & Maschio (1971) and Bonazzola & 
Schneider (1974), Bonazzola et al. (1993) present a formulation which allows for the most general 
stress-energy tensor consistent with a spacetime having the properties of stationarity, axisymmetry, 
and circularity. The metric for such a spacetime can be expressed as 



2. 



FORMALISM 




(1) 



gapdx'^dxl^ = -e^^dt"^ + e'^^G^^ sin^ e{d(l) - N'^dtf + e^^'^-''\dr'^ + r'^dO'^), 



(2) 
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where the metric functions u, G, , and C are functions of (r, ^). (A spacetime having the prop- 
erties of stationarity and axisymmetry, but not circularity, would have an additional off-diagonal 
term in the metric (Bocquet et al. 1995; Carter 1973)). For the spacetime to have the property of 
circularity in addition to the properties of stationarity and axisymmetry, it is necessary that the 
electromagnetic current 4-vector and fluid 4-velocity be parallel to a general linear combination 
of the Killing vectors (Carter 1973). For example, in the coordinates of equation (2), J* and J"^ 
would be the only nonvanishing components of the electromagnetic current 4-vector. A further 
consequence of this is that At and Aff, are the only nonvanishing components of the electromagnetic 
potential 1-form (Carter 1973). 

From the Einstein equations, Bonazzola et al. (1993) derive a Poisson equation for each of the 
metric variables, with source terms that depend on the metric variables and on the components of 
the stress-energy tensor: 

(3) 
(4) 
(5) 
(6) 

where 

= rsinON'^, (7) 
G = rsinOG, (8) 

and A2, A3, and A3 are respectively the 2D flat space Laplacian, the 3D fiat space Laplacian, and 
the (f) component of the 3D flat space vector Laplacian: 

A = — + + (9) 

^ dr"^ r dr dO'^ tan 6 89 ' 

A3 = A3- ^ ^ (11) 
sm^ 9 

The source terms in equations (3-6) involve the metric functions and components of the stress- 
energy tensor: 



A3Z/ = 




A3iV'^ = 




A2G = 




A2C = 





AttGn e^^^-"^ {E + S\) + ]- e-^^GV2 sin^ 9 {dN't'f - du 9(ln G), (12) 



= - ^^'"^^f^^ ^.\ -rsin9 dNt> d [in {e-'^G^) ] , 



Q2 J, gjj^ Q 



(13) 



(7/^ = 87rG7ve2K-'^)Grsine(5''^ + S'^e), (14) 



(7^ = 87rGive2«-'^)5*^ + ^e-^"GVsin2 0(aAr^)2-(5z/)2. (15) 
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In these expressions the notation dX dY denotes 



dXdY IdXdY 



In addition, Gat is Newton's constant, and the contributions from the stress-energy tensor are^ 

E = T^fsn^n'^, (17) 
la = -KpUpT^P, (18) 

Safi = Kphfi^TP'', (19) 

where n° is the unit four vector orthogonal to the spacelike time slices and hap = QajB + nanp is 
the spacelikc projection tensor. In the present coordinates tt-q, = {—N, 0, 0, 0), where iV = e*^ is the 
lapse function. For the stress-energy tensor of a perfect fluid [the first part of equation (1)], 

E^^ = r\e + p)-p, (20) 
(JP^)^ = e-'^G r sine {E^^ +p)U, (21) 

(5^^); = p, {S^n\=P^ {S^^)\ = p+{E^^+p)U\ (22) 

where T = (1 — C/^)~^/^ and U = e~'^'^Gr sin 6{Q-N'f'), and the superscript "PF" stands for "perfect 
fluid." In this expression for U, the scalar ft is the angular velocity as measured by a static observer 
at infinity; it relates the nonvanishing components of n° through the equation u'^ = Q.u^ For 
the standard stress-energy tensor of the electromagnetic field [the last two terms of equation (1)], 
restricted to a poloidal field, we have 

^EM ^ ^{EiE' + BiB'), (23) 

(jEM)^ ^ l.e^'^-^^Gr'^smeiE'B^-E^B'), (24) 

(.S™); = ^{EeE<'-ErE' + B0B'-BrB^), (25) 

in which the superscript "EM" identifies the electromagnetic contributions. In these expressions Ei 
and Bi are the components of the electric and magnetic fields as measured by an Eulerian observer 
(i.e. an observer with four velocity n°): 

Ea = n^Fap (27) 



^We do not use E to denote the magnitude of the electric field, although we will use the subscripted notation Ei 
to refer to components of the electric field. 

^As will be discussed below, in the limit of infinite conductivity ("frozen-in magnetic fields"), stationary stars with 
magnetic fields must be uniformly rotating. 
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dAt ..^dA, 



dr 



dr 



dAt ,.^9 A, 



= 0, 



2 ^a/Spa"^^ 

e" OA, 



OA, 



Gr^ sine 86 
where ea/3pa is the Levi-Civita tensor. 



Gsin^ dr 



de 



-,o 



+ AT'/'. 



de 



,0 



(28) 
(29) 
(30) 



The quantities At and ^0 are also determined by Poisson equations, which derive from the 
Maxwell equations in curved spacetime (Bocquet et al. 1995): 



AsAt 



(31) 
(32) 



where A'^' = A^p/ {r sin 9) . The sources can be expressed 



a A, = -47r e2(^-'^) (^^ J* + J-^) + e'^'^gt^ dAt dN't' _ (2 + g-^-^,,) dA^dN'l' 



-idAt + 2iV^a^,) a [m (e-G)] - ^ + 
-47r e^^-^'^GVsine (^J"^ - Ar<^J*) + e-^'^GVsin^SAT'^ (dAt + N^dA^ 



dA, 



a 



A'l' 



rsva.6 



dA^d [ln(e-2^G)] . 



(33) 



(34) 



A theorem of Cowling (1934) states that an axisymmetric magnetic field cannot be generated 
or maintained by the motion of a fluid. The theorem relies on the fact that finite resistivity involves 
dissipation, leading to magnetic field decay.^ Hence stationary models of neutron stars in magnetic 
fields require a separation of dynamical and dissipative time scales, encoded in an assumption of 
infinite conductivity [magnetic fields "frozen in" and carried with the fluid, a common assumption 
in astrophysics (Alven 1950)]. This assumption is exceedingly well justified for neutron star matter, 
the Ohmic dissipation time scale being larger than the age of the universe (Goldreich k, Reisenegger 
1992). Cowling's theorem is thus effectively nullified. In addition, as shown by Bonazzola et al. 
(1993), the assumption of infinite conductivity leads to the requirement of uniform rotation, as well 
as the relation 



At = —^Aff, + const 

inside the star, where the constant is determined by the total electric charge of the star. 



(35) 



Closure of the system of equations requires relations involving some quantities appearing in 
the source equations (12-15) and (33-34): the rest-frame energy density e and pressure p, and 



''Ambipolar diffusion and Hall drift (Goldreich & Reisenegger 1992) will actually dominate magnetic field decay 
in magnctars, but even for the ultra strong fields considered here the decay time will be hundreds if not thousands 
of years. 
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the components J* and J'^ of the electromagnetic 4-current. For a uniformly rotating stationary 

;/3 



star in a magnetic field, local energy-momentum conservation {T°'^.q = 0) yields the equations of 



stationary equilibrium: (Bonazzola et al. 1993) 

-^p,i + y,i-{\nT),,--^Fi^r = Q, (36) 
e + p e + p 

in which = dX/dx^. Derivatives with respect to the coordinates = {r,6) give the only 
nontrivial equations, due to the symmetries of the problem. We recall that v = (1/2) \n{—gtt) here 
plays a role like that of the gravitational potential in the Newtonian case, and that T is the Lorentz 
factor associated with the fluid's rotational velocity. By analogy with the Newtonian case, the 
terms in equation (36) may be thought of as (from left to right) the pressure force, gravitational 
force, centrifugal force, and Lorentz force. 

For a one-parameter equation of state, e = e(n), p = p{n), there is a first integral of the first 
term in equation (36): 

^(^) = r r /^ ! r A 7^ dn', (37) 
Jq e[n') + p[n') an' 

where we will assume that h{0) = corresponds to the surface of the star.^ This, together with 
the adoption of a "current function" f{A^), 

-^{jt>-nj')=f{A^), (38) 
e + p 

give the equations of hydrostatic equilibrium a first integral (Bonazzola et al. 1993) : 

h{r, 9) + zy(r, 6) - In r(r, 9) + M(r, 9) = C = constant, (39) 
where the "magnetic potential" M(r, 9) is given by 

M{r,9) = M{A^{r,9)) = - / f{x)dx. (40) 

Jo 

While one is free to choose the current function f{Aff,), equation (38) represents a significant restric- 
tion on the form of the electromagnetic current that allows the existence of stationary solutions. 
The constant C is determined by an input parameter, e.g. the density specified at some point in 
the star. 

In summary, the formalism of stationary neutron stars in poloidal magnetic fields with a one- 
parameter equation of state consists of a closed system of eleven variables [four metric variables, 
energy density, pressure, two components of the electromagnetic potential, two components of the 



®For example, at zero temperature and in chemical equilibrium, the pressure and energy density are functions 
only of the baryon density n. By virtue of the first law of thermodynamics, h{n) is seen to be the logarithm of the 
enthalpy per baryon: h = ln[(e + p)/n] + constant. 
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electromagnetic current, and the "heat function" h of equation (37)]; eleven equations [four Poisson 
equations (3-6) for the metric variables, two Poisson equations (31-32) for the components of the 
electromagnetic potential, the relation (35) between the components of the electromagnetic poten- 
tial, the equation of state, the relation (37) between the heat function and e and p, the first integral 
(39) of the equations of hydrostatic equilibrium, and the restriction (38) on the electromagnetic 
current]; three input parameters (angular velocity, total electric charge, and maximum density); 
and one input function [the /(^</,) in equation (38)]. 

3. WHEN IS THE MAGNETIC FIELD TOO LARGE? 

As noted in the introduction, in this study we restrict ourselves to the static case. This involves 
a number of simplifications, including the vanishing of N^^, At and J*, and the absence of surface 
charges generally present on perfectly conducting rotating bodies. 

Bocquet et al. (1995) present numerical calculations aimed at determining "the maximum 
mass configuration among all static magnetized models" for several equations of state and for the 
choice of constant current function /(A^) = /q. These static configurations are determined by two 
parameters, which Bocquet et al. (1995) took to be the central density and /q. They considered 
sequences of constant magnetic dipole moment M, which is defined in terms of the asymptotic 
behavior of (the orthonormal components of) the magnetic field, 

{2M cos 9/r^) = B^^) \r^oo = {e^'^-'^/Gr^ sin 9) (OA^/ de)\r^oo , (41) 

and then determined the maximum mass for each value of M. For example, using the "Pol2" 
7 = 2 polytropic equation of state (Salgado et al. 1994), they reported the maximum mass among 
all static configurations to be 4.O62M0, with a magnetic moment of 1.122 x 10^^ A m^. We have 
calculated, also using the Pol2 EOS, a similar configuration, which is pictured in Figure 1, and 
may be compared with Figures 5 and 6 of Bocquet et al. (1995). We point out that in Figures 1 
through 6 that we use the same polytropic constant in the Pol2 EOS as that employed by Bocquet 
et al. (1995). However, in later figures and in all tables, we have altered the polytropic constant to 
a more realistic value so that the maximum mass of static configurations is 2.OM0, and then scaled 
the results from the other groups accordingly in these tables and figures. 

Bocquet et al. (1995) also note that "For magnetic fields higher than [this configuration], no 
stationary configuration can exist and the numerical procedure. . . fails to converge." Note that 
this is not a question of stability — nothing is being claimed about the stability of the high field 
configurations they achieve — the question pertains to even the existence of stationary solutions. 
The fact that the maximum mass and magnetic moment are determined by the failure of their code 
to converge leaves one wondering whether a genuine physical limit has been reached, or whether 
the result simply represents the failure of the numerical method to find solutions that do in fact 
exist. 
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This situation is in contrast to the case of rotation: It is well known that there is an upper 
limit to the angular velocity of uniformly rotating stars. Numerical codes that attempt to solve for 
stellar configurations with uniform angular velocities above this limit will fail to converge. However, 
there is also a well-defined physical basis for this "Kcplcrian" limit, namely, mass shedding.^ 

At the mass shedding limit, a fluid element at the equatorial surface undergoes geodesic motion: 
it can remain in that orbit without any pressure support from the star. Thus, it is not necessary 
to rely only on the failure of the code to converge to determine the mass shedding limit; one can 
quantitatively test for geodesic motion at the equatorial surface. 

The approach to the mass shedding limit can also be visualized using the first integral of 
the equations of hydrostatic equilibrium — the relativistic generalization of the Bernoulli equation. 
Prom equation (39) , the "Bernoulli equation" for a uniformly rotating star without a magnetic field 
is 



Figure 2 demonstrates how the Bernoulli equation can be used to visualize the approach to the mass 
shedding limit. In this figure the dashed lines represent the dot-dashed lines represent — InF, 
and the thick solid lines represent v — InF, all as functions of a compactified radial coordinate in 
the equatorial plane. The dotted lines represent C, so that h is given by the distance between the 
dotted and thick solid lines. The surface of the star {h = 0) then corresponds to the intersection of 
the dotted and thick solid lines. In the lower panel of Figure 2, which represents the mass shedding 
limit, the equatorial surface of the star coincides with the maximum of — InF. For larger angular 
velocities, C would be larger than this maximum, and there would be no surface of the star at finite 
radius. 

Similar plots can be made for stationary stars with magnetic fields. For a nonrotating star in 
a poloidal magnetic field, equation (39) reduces to 



Now, in the equatorial plane, along which direction are magnetic forces exerted? The answer 
depends on the direction of the magnetic field, which in turn depends on the current distribution 
in the star. From equation (38) (with $7 and J* set to zero in the present static case), the electric 
current density is proportional to {e+ p), so that in the typical case one would expect the current 
to vanish on the surface of the star. The current measured by a local observer must also vanish on 
the axis of symmetry. Hence the current J^^^ measured by a local observer in the equatorial plane 
is expected to have a structure that peaks somewhere inside the star and vanishes at the origin and 



In reality, the rotation may be more severely limited by a gravitational instability to non-axisymmetric pertur- 
bations. The instability might, however, be damped out by bulk and shear viscous effects (Lindblom, & Detweiler 
(1977); Imamura, Durisen, & Friedman (1985); Friedman, Ipser, & Paxker (1986); Ipser, & Lindblom (1989); Sawyer 
(1989)) To the extent that the Keplerian limit represents a reasonable estimate of the upper limit of rotation, our 
analysis here illuminates the stark contrasts that exist between the instabilites caused by rotation and magnetic fields. 



h{r, e) + v{r, e) - lnr(r, 6) = C = constant. 



(42) 



h{r, e) + u{r, e) + M{r, 6) = C = constant. 



(43) 
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the surface, as shown in the left panel of Figure 3 for the Pol2 EOS. Since the magnetic field lines 
tend to circle around a point in the vicinity of the maximum current (see the right panel in Figure 
1), in the equatorial plane the field reverses direction inside the star. Accordingly, the Lorentz force 
reverses direction inside the star, as displayed in the right panel of Figure 3. 

The Lorentz force and the other "forces" acting in the equatorial plane are derivatives of the 
quantities appearing in equation (43). These quantities are plotted in Figure 4 for the Pol2 EOS, 
which is similar to Figure 2, except that the dot-dashed lines now represent the magnetic potential 
M{r,9) in equation (39), given in this case by equation (40) with / = constant = /q. In the inner 
portion of the star, the Lorentz force behaves like a centrifugal force, pushing outward and allowing 
the star to support more mass. As seen in the lower panel of Figure 4 (and in the left panel of 
Figure 1), this outward force can be strong enough to cause the maximum of h (and hence e) to be 
off-center. However, the reversal of the force at larger radii makes the total "potential" confining. 
Thus, in contrast to the case of rotation, there is no mass shedding limit at the equatorial surface 
that clearly indicates the nonexistence of stationary solutions and explains the failure of the code 
to converge for large magnetic fields. 

If mass shedding does not occur, is there some other identifiable physical cause that prevents 
the existence of stationary solutions for sufficiently large magnetic fields? Bocquet et al. (1995) 
note that for sufficiently large fields, the total (fluid + magnetic) stress tensor has a component 
on the symmetry axis that goes from being positive (pressure) to negative (tension), causing the 
star to have a characteristic "pinched" shape (see the left panel of Figure 1). This occurs when the 
"magnetic pressure" exceeds the fluid pressure. They note that in the largest mass configurations 
they obtain, the ratio of the magnetic to fluid pressures at the center of the star is of order unity. 
However, there are two reasons that argue against considering this as a definitive physical reason 
preventing the existence of stationary solutions. The first is a matter of principle: the equilibrium 
of fluid elements relies on a balance between gravity and and the gradients of stresses; the absolute 
sign of the stresses themselves is not of fundamental importance. The second reason is the results 
observed in practice: the value of the ratio of magnetic pressure to fluid pressure at the center of 
the star reported by Bocquet et al. (1995), for the putative maximum mass conflgurations, varies 
quite noticeably among different equations of state, and its value is not predictable. This is not 
what one would expect of a precise physical criterion for the nonexistence of stationary solutions. 

For Newtonian stars, Chandrasekhar &; Fermi (1953) identifled gravitational binding, or a 

negative total energy (excluding rest mass), as a necessary criterion for the dynamical stability of 
equilibrium solutions. For nonrotating polytropes, they used a generalization of the virial theorem 
that accounts for magnetic fields to show that a negative total energy requires that the magnitude 
of the gravitational potential energy must be greater than the magnetic field energy. 

One might wonder whether a criterion of gravitational binding could represent a physical 
upper limit on stellar magnetic fields in the relativistic case as well. While total gravitational 
and magnetic energies are not well defined in relativity, in some sense the gravitational mass M 
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contains all forms of "energy," including "magnetic energy." Since the total baryon number [or 
total baryon mass Mb = (baryon mass) x (baryon number)] is a well defined quantity, one might 
consider using the relation M — Mb < as a criterion for gravitational binding. As it turns 
out, the putative maximum mass and maximum field configurations reported by Bocquet et al. 
(1995) are still "gravitationally bound" by this criterion. The question remains: As the influence 
of the magnetic field increases, is this apparent transition from existence to nonexistence of static 
solutions, occuring at finite "gravitational binding," a physical result or a numerical artifact? 

In order to explore these questions we have searched for the maximum mass in a different way 
than Bocquet et al. (1995). As mentioned previously, Bocquet et al. (1995) computed sequences of 
constant magnetic moment Ai by suitably adjusting the central log-enthalpy he and the value /o of 
the (constant) current function /. For each value of Ai they determined the maximum mass. The 
overall maximum mass was obtained with the largest value of A4 for which convergence could be 
achieved. In contrast, we have chosen a more direct means of exploring parameter space. For each 
value of maximum log-enthalpy /tmax (which for large fields will not coincide with he), we found 
the largest values of /o for which the code converged. Specifying /imax instead of he allows for the 
possibility of vanishing density at the origin, i.e. toroidal configurations. 

A close inspection of the forces in configurations near the failure to converge reveals an apparent 
physical cause, for the choice of constant current function, of the failure to find stationary solutions 
for sufficiently large magnetic fields: When a sufficient quantity of matter has been pushed off-center 
by magnetic forces that gravitational forces begin to point radially outward in the equatorial plane, 
a topological change to a toroidal configuration ensues. This is illustrated by Figure 5 for the Pol2 
EOS, which shows the late stages of iteration of a configuration with a value of /o that is slightly too 
large for convergence. We emphasize that the stages depicted in this figure are not valid stationary 
solutions to the Einstein equations; neither do they represent a true evolution. Nevertheless, the 
sequence is suggestive of possible dynamical outcomes: a transition to a toroidal topology, expansion 
of the torus to large radii, and increasing compactification of the toroidal configuration of matter. 
As the iterations proceed, the outward pointing gravitational force (positive —du/ds or —dv/dr) 
becomes more and more pronounced. The central evacuation begins, however, with a tiny outward 
gravitational force at some radius. 

Since we were unable to find any convergent toroidal solutions for the case of a constant 
current function /, we can therefore identify a quantitative criterion for the boundary of existence 
of stationary solutions in this case: this boundary is characterized by du/dr = at some off-center 
location in the equatorial plane. By symmetry, it is always the case that du/dr = at the origin; 
furthermore, the positive-semidefiniteness of energy density and pressure ensure that d'^vjdr'^ > 
at the origin (this can be verified from the Green's function expansions of the metric functions 
presented in Appendix A). This means that the critical condition du/dr < cannot occur at the 
origin or an infinitesimal region surrounding it, but our numerical calculations indicate that it first 
occurs very close to the origin. 
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These matters are illustrated in Figure 6 with the Pol2 EOS, which shows the radial profile 
of the gravitational force for increasing values of /o — at a particular value of maximum density — 
including the largest value of /o for which convergence was achieved. In the upper panel, a flattening 
of the gravitational force near the center with increasing magnetic field is apparent. The insets in 
the upper panel show the region near the origin. The lower left inset shows —dv/dr as computed 
with a formula obtained by differentiating equation (A2). As required by the condition d'^vjdr^ > 
at the origin, according to this formula the first off-origin grid point has negative value of —dv/dr. 
For the largest value of /o, the next grid point is less negative, but still does not reach our asserted 
condition —du/dr = 0. However, the upper right inset in the upper panel shows values of —du/dr 
near the origin computed from a the centered finite difference formula, which is closer to what the 
discretized configuration actually "feels." Computed in this way, the first off-origin grid point has a 
barely positive value of —du/dr for the largest value of /o for which convergence was achieved. The 
lower panel of Figure 6 shows the gravitational force for the largest value of /o, for three different 
grid resolutions. The insets of this panel again compare the "analytic" and "numerical" derivatives. 
These computations show that with increasing resolution, the "analytic" derivative near the origin 
gets closer and closer to the critical condition —dv/dr = 0. In practice, however, the maximum 
mass and magnetic field configuration for a given stellar maximum density can be identified by the 
appearance of a small positive value of —dv/dr, as computed with a centered finite difference. 

For the values of maximum density we studied, and for a constant current function /, we 
did not find any stationary toroidal solutions. The toroidal configurations continued to expand in 
radius and compress into thinner and thinner rings until the region covered by matter consisted 
of only a few gridpoints, at which point the code would fail. A determination of the outcome of 
the evolution of such configurations — whether to a stationary solution characterized by a different 
current function, dispersal to infinity, or even the formation of a toroidal event horizon (Shapiro, 
Teukolsky, &; Winicour 1995) — would appear to require a fully relativistic evolution code. 

It remains to be seen if these results — namely, the prescription for determining when magnetic 
forces are too strong for the existence of static configurations, and the lack of converged toroidal 
solutions — will hold for more general current functions. 

4. SEQUENCES OF CONSTANT BARYON MASS AND CONSTANT 

MAGNETIC MOMENT 

Our computations of static neutron stars with Tiltra-strong magnetic fields determined by a 
constant current function are summarized in Table 2 and Figures 7 and 8. Figure 7 contains some 
of the older EOSs employed by Bocquet ct al. (1995). These include BJI, which is model IH 
of Bethe &; Johnson (1974) [we derived our table from that listed in Malone, Johnson, & Bethe 
(1975)], and PandN, which is from Pandharipande (1971). Results for more recent EOSs are 
displayed in Figure 8. These include Akmal, which is from Akmal et al. (1998) and is based on 
a potential model description of dense matter and represents the most complete study to date in 
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which many-body and special relativistic corrections are incorporated. PCL is taken from Prakash 
et al. (1995), and is based on a relativistic field-theoretical description of dense matter starting 
from the Lagrangian proposed by Zimanyi & Moszkowski (1990). This approach easily allows for 
the inclusion of additional softening components: the case in which hypcrons are present is labelled 
PCLhyp. For all EOSs except Pol2, we employed the EOS of Baym, Bethe, k Pethick (1971) and 
Baym, Pethick, & Sutherland (1971) at densities below about 1/2 the nuclear saturation density. 

The rationale for exploring a wide variety of EOSs, even some that are relatively outdated, is 
two-fold. First, it provides contrasts among widely different theoretical paradigms. Second, it illu- 
minates general relationships that exist between the pressure-density relation and the macroscopic 
properties of the star such as the maximum mass and the radius. 

In these figures the lower thick solid curves show the gravitational mass M as a function of 
radius for static stars without magnetic fields, which are spherical. The upper thick solid curves 
show the outer boundaries beyond which no static solutions were found (sec §3). The lighter solid 
curves are sequences of constant baryon mass Mb (and varying magnetic moment Ai), while the 
dotted curves are sequences of constant magnetic moment (and varying baryon mass). The lighter 
shaded regions indicate configurations in which the magnetic field is sufficiently strong that the 
maximum density is off-center. The small slivers of darkly shaded regions towards the left sides of 
the plots indicate solutions which are gravitationally unbound (M — Mb > 0) ; these are expected 
to be dynamically unstable. 

As with rotation, magnetic fields allow neutron stars with a particular EOS and baryon number 
to have larger masses and equatorial radii compared to the field-free case. In Figures 7 and 8, the 
configurations of maximum mass that can be reached by uniform rotation (without magnetic fields) 
are marked with an "X" . For all the EOSs displayed, our result for the maximum mass attainable 
with a magnetic field governed by a constant current function is noticeably larger than that attained 
by rotation. This contrasts with the results reported by Bocquet et al. (1995) for the maximum 
mass attainable with magnetic fields, shown in Figure 7 with crosses. Out of the three EOSs we have 
in common, only in the case of the polytropic EOS do Bocquet et al. (1995) obtain a significantly 
larger mass with magnetic fields than with rotation. 

In the absence of significant accretion, constant baryon mass sequences are of interest as po- 
tential evolutionary paths. This is well motivated in the case of uniformly rotating non-magnetic 
stars: As angular momentum is slowly dissipated by gravitational radiation, the star moves along 
a sequence of constant baryon mass until it either stops rotating (for "normal" sequences termi- 
nating on the spherical mass vs. radius curve) or collapses to a black hole [for "supramassive" 
sequences which exist solely by virtue of rotation; stars on such sequences may exhibit the interest- 
ing phenomenon of spin-up during angular momentum loss just before collapse to black hole (Cook, 
Shapiro, & Teukolsky 1992, 1994a,b; Salgado et al. 1994)]. These scenarios are well motivated be- 
cause of the expectation that once viscosity brings a dynamically stable star into uniform rotation, 
it will not spontaneously begin to differentially rotate as it loses angular momentum. 
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On the other hand, representing evolutionary sequences by the constant baryon mass sequences 
for magnetized stars pictured in Figures 7 and 8 is an oversimpHfication. Even though a star's mag- 
netic field will slowly (on dynamical time scales) decay due to Hall drift and ambipolar diffusion, 
there is no guarantee that the star's configuration will proceed along the paths pictured in the 
figures. This is because the case of magnetic fields is more analogous to differential rotation rather 
than uniform rotation, the necessary choice of a current function in the magnetic case [see equa- 
tion (38) and surrounding discussion] corresponding to the choice of a rotation law in the case of 
differential rotation. As a star's magnetic field decays, it is not obvious that its current function 
will remain the same.^ Perhaps the study of several different current functions could shed light 
on probable evolutionary sequences. For example, an analysis of how the mass varies with the 
functional form of the current function (at fixed baryon mass and magnetic moment) could give 
an idea of how the slow evolution with magnetic field decay might proceed. There is of course no 
principle of "conservation of magnetic moment," but since the time scale for magnetic field decay 
is slow, this procedure seems like a plausible opening exploration. 

In connection with the constant baryon mass sequences, we here comment on a curious feature 
in Figures 7 and 8. For the potential model EOSs BJI, PandN and Akmal, the topology of these 
sequences near the maximum spherically symmetric star appears to be different than for the poly- 
tropic or relativistic field theoretical models PCL and PCLhyp. In the former, there is a minimum 
in the baryon mass above the spherical, non-magnetized, sequence. To determine if this feature 
was related to the possible acausal behavior of potential models at high density, we modified the 
PandN EOS to go over to a causality limit EOS when necessary, but found that the topology was 
unaltered. Instead, the effect appears to be related to the fact that all forms of energy contribute to 
the magnetic field: While the Newtonian intuition (and the relativstic behavior at lower densities) 
is that magnetic fields always increase the gravitational mass of a star of given baryon number, the 
fact that the energy density in magnetic fields can be a nontrivial source of gravitation means that 
this self-gravitating tendency can compete with the tendency of the Lorenz force to help support 
the star. Perhaps this occurs near the spherical maximum mass in the case of the potential model 
EOSs, with the result that configurations near the spherical maximum mass have nowhere to go 
but down in gravitational mass when magnetic fields (governed by a constant current function) 
are added. If a star did preserve its current function and follow a supramassive constant baryon 
mass sequence possessing such a gap, upon reaching the minimum mass of the sequence it could 
catastrophically collapse to a spherical neutron star rather than a black hole. While the possibility 
is admittedly remote, this could be a novel form of energy release in a baryon- free environment, 
giving rise to a mini-7-ray burst. 

Except for the spherical stars, those with no magnetization indicated by the lower solid line 
in these figures, the stability of the configurations has not been studied. For spherical stars, it 



'^Note that constant current function is not analogous to uniform rotation. Formally, the analogue to uniform 
rotation would be uniform magnetic vector potential, which of course means no magnetic field. 
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is well-known that the configurations with larger radii than the maximum mass configuration are 

stable. One might speculate that constant baryon mass sequences which terminate on the stable 
side of the spherical AI vs. R curve are stable, while those ending on the unstable side would be 
unstable, and that for supramassive constant baryon mass sequences the minimum mass marks 
a change from stability to instability. But this remains to be determined. As with differentially 
rotating stars, it is necessary to do a normal mode analysis, or even a fully relativistic evolution; see 
Baumgarte, Shapiro, & Shibata (2000). Incidentally, it is interesting to note that these authors find 
a differentially rotating configuration with a "red blood cell" shape similar to our extreme magnetic 
configurations, and that this configuration is dynamically stable. Of course, their configuration was 
not subject to MHD instabilities that may come into play (Spruit 1999a,b). 

5. SUMMARY AND OUTLOOK 

In summary, we present a method of computing the structure of axisymmteric relativistic stars 
that combines elements of previous approaches, and report tests of our code. A quantitative method 
of determining the outer envelope (in the mass vs. radius plane) of configurations attainable with 
poloidal magnetic fields governed by a constant "current function" [see equation (38)] has been 
found: magnetic fields are too large for static configurations to exist when the magnetic force 
pTishcs a sufficient amount of mass off-center that the gravitational force points outward near the 
origin in the equatorial plane. (In our coordinates an outward gravitational force corresponds 
to —du/dr > 0.) We obtain larger masses of neutron stars in ultra-strong magnetic fields than 
have been reported previously for various equations of state, and performed computations with 
three representative modern EOSs. Sequences of constant baryon mass and constant magnetic 
moment are displayed. For all EOSs studied, the maximum attainable mass of static stars with a 
magnetic field determined by a constant current function is noticeably larger than that attainable 
with uniform rotation and no magnetic field. 

The results presented here are only an initial step in exploring possible configurations of neu- 
tron stars with strong magnetic fields. As we mention below, configurations with azimuthal field 
components will be of physical interest, which implies that three-dimensional geometries should 
also be considered. But even with attention restricted to poloidal fields, we have only scratched the 
surface of possible configurations, as we have only considered a single current function. Bocquet 
et al. (1995) make brief mention of computations with a few other current functions. We have 
performed a handful of exploratory computations using a polytropic EOS and other current func- 
tions and have found some toroidal solutions. These toroidal solutions were not attainable with the 
computational approach of Bocquet et al. (1995), since their method involved the specification of a 
finite density at the center of the star. In the case of toroidal configurations, the simple condition 
determining the boundary of existence of static configurations will have to be generalized, since 
there is no matter at the center. These explorations will be reported in detail elsewhere. 

Our work here has focused on the effects magnetic fields have on general relativistic structure. 
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ignoring the effects of intense magnetic fields on the EOSs. Recently, the direct effects of magnetic 
fields on the EOS have also been investigated (Chakrabarty 1996; Chakrabarty, Bandyopadhyay, 
&: Pal 1997; Yuan & Zhang 1999; Broderick, Prakash, & Lattimer 2000). Substantial effects on the 
EOS above nuclear saturation densities are generally produced by fields in excess of 10^^ G, which 
is of the order of the maximum central field strengths found in this paper. The generic effects 
on the EOS include softening due to Landau quantization, which is, however, overwhelmed by 
stiffening due to the incorporation of the magnetic moments of the various particles in neutron star 
matter (Broderick, Prakash, & Lattimer 2000). (Note that the important S^/Svr term is already 
included in our study.) Work is in progress (Cardall et al. 2000, in preparation) to provide fully 
self-consistent calculations of neutron star structure, in which the direct effects of magnetic fields 
on the EOS will be included in addition to the structural effects considered in this work. 

An issue which we defer to future work is the question of stability. For systems governed 
by a finite number of parameters, a generalization (Sorkin 1982) of the familiar one-dimensional 
turning point method can be employed; see Friedman, Ipser, &; Sorkin (1988) for an application to 
uniformly rotating relativistic stars. However, as with differentially rotating stars, this generalized 
turning point method is not really applicable in the present case. This is because the need to 
specify a current function (or rotation law in the case of differential rotation) means that defining 
a particular configuration requires the specification of an infinite number of parameters. 

Another issue that needs further explication before the physical relevance of the results pre- 
sented here can be fully assessed is the generation of magnetic fields. We have mentioned the 
mechanism of Duncan &; Thompson (1992), the generation of magnetic fields during the smoothing 
of differential rotation. However, the azimuthal dragging of field lines by differential rotation leads 
to nonvanishing azimuthal field components, in constrast to the poloidal fields studied here. It is 
not clear whether fields with azimuthal components would evolve into poloidal configurations, or 
whether there are mechanisms to directly generate strong poloidal fields. It would be of interest 
to explore the possibility of finding stationary solutions with toroidal magnetic fields, and in three 
dimensions. While this would involve more nonzero metric components, perhaps methods similar 
to those employed in this paper could be employed; see Bonazzola, Gourgoulhon, & Marck (1998). 

We wish to thank E. Gourgoulhon for helpful communications concerning the calculations 
of Bocquet et al. (1995). We are grateful to Dany Page and Ralph Wijers for their help in the 
preparation of Table 1. Research support from DOE grants FG02-87ER40317 (for CYC and JML) 
and FG02-88ER-40388 (for MP) are gratefully acknowledged. 

A. NUMERICAL PROCEDURES AND TESTS OF THE CODE 

Bonazzola et al. (1993) have developed a "spectral method" to solve equations (3-6). This 
method involves expanding the solution on a set of basis functions having the analytical properties of 
the solution. They also try to choose basis functions for which there exist fast transform algorithms. 
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Their method involves the use of two grids, an "inner" region including the origin and an "outer" 

region that reaches to infinity. Instead of this "spectral method," we have chosen to solve the Poisson 
equations using Green's functions, similar to the methods of Komatsu, Eriguchi, &: Hachisu (1989) 
and Cook, Shapiro, k. Teukolsky (1992, 1994a,b). It is convenient to compactify the radial domain 
0<r<ootoO<s<l via the change of variables 

. = ij(^), (Al) 

where R is some length scale. Taking into account the azimuthal and equatorial symmetries, and 
imposing the boundary conditions [all metric functions finite at the origin; {y, N'^', C)|r^oo ~^ 0, 
G|r-+oo — ^ 1], we find equations (3-6) to yield 
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in which the quantities on the right hand side are given by equations (12-14). The symbols Pn{x) 



-19- 



and P^{x) denote the Legendre polynomial and the associated Legendre function, respectively. 
The source ct^ requires special consideration. The Green's function of the 2D Laplacian has a 
Inr term, as is apparent from equation (A5). This term must vanish as r ^ cx) in order to have 
vanishing boundary conditions at infinity. This is not a problem for equation (5) for G, since the 
sm.9 factor in the source Gq [see equation (14)] guarantees that the Inr term vanishes everywhere. 
However, there is no such factor in the source cr^ [see equation (6)] . A genuine solution to equation 
(6) satisifying the boundary condition Clr^oo will have a hir term that vanishes as r — > oo; 
hwt in the intermediate steps of an iterative procedure to solve the nonlinear equations there is no 
guarantee that this will be so, leading to a potential instability. 

Bonazzola et al. (1993) have a resolution of this difficulty which we adopt here. In order for 
the Inr term to vanish as r — oo it is necessary that 

/ / a^{r,e)rdrde = Q. (A9) 



This condition is called the "virial theorem" by Bonazzola et al. (1993). In terms of the variable s 
it can be written as 

^^c7c,o(5) = 0, (AlO) 



(1 - Sf 

where 

/•27r 

ac,o(s)= / dOa^is^O). (All) 

The trick is to divide the source into two pieces. One piece, u™, contains the "matter terms" 
(those involving components of the stress-energy tensor); the other piece, a^, contains the "field 
terms," those involving only the metric variables. The virial theorem can then be written 



f'^ dss ^ dss f 



(A12) 



This equation will be satisfied for the actual solution to the Einstein equations, but will not be 
satisfied in the intermediate steps of the iteration procedure. To avoid the potential logarithmic 
divergence associated with the failure to satisfy equation (A12), the source cr^ is replaced by cr"* + 
Xa^, where 



A = - 



/ 



ds s f , ^ 

^r37)3-c,o(-) 
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In this way equation (AlO) is satisfied at each step of the iteration — avoiding the potential loga- 
rithmic singularity — but with A / 1 in the intermediate steps. At the end of the iteration process 
A must approach 1 for the computed metric functions to represent a valid solution to the Einstein 
equations. Finally, in equation (A5) is given hy = Rr {a^ + Acr^). 

It is convenient to know ahead of time the location of the equatorial surface on the grid. To 
achieve this we employ a scheme like that Bonazzola et al. (1993) use to divide their computational 
domain into "inner" and "outer" grids. Specifically, we specify that the equatorial surface be 
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located at the radial position s = 0.5. Prom equation (Al), this makes the equatorial radius equal 

to R. Since R is some chosen constant, this involves a nonstandard system of units. Operationally, 
the value of Newton's constant Gjsf is adjusted at each iteration in such a way that s = 0.5 docs 
indeed correspond to the equatorial radius of the neutron star, identified by /i = 0, where h is 
defined by equation (37). The physical value of the equatorial radius scales as ^/G^; other physical 
quantities also involve various powers of Gjv- 

The scheme is implemented as follows. As with cr^ described above, the source term is 
divided into two parts, with the "matter part" crJJ* containing source terms deriving from the 
stress-energy tensor, and the "field part" at containing terms involving derivatives of the metric 
variables. The Poisson equation for u is solved in two parts: Aai/"* = a'^, where a'^ = a"^ jG^\ 
and ^j,v^ = at- The full value of v is then v = G^i'™' + z^-^. The demand that h vanish at 
s = s* = 0.5 in the equatorial plane, together with equation (39), yields the appropriate value of 
Newton's constant at each iteration: 

^ _ + - InP + M) U=s^^^ - (^/ - InP + M) 

-.ml --"1.1 ■ V / 



— j/i' 

IS — 5* ^ \S — 5max 

In this expression, h\s=s^^^ is an input parameter that, via equation (37) and the equation of state, 
specifies the maximum density in the equatorial plane. This is not necessarily the central density; 
its location Smax must be determined at each iteration. Specifying the maximum density while 
allowing its location to "float" allows for the possibility of toroidal configurations, a possibility 
excluded by the method of Bonazzola et al. (1993) and Bocquet et al. (1995). 

In order to test the ability of our code to solve the Einstein equations for axisymmetric configu- 
rations, we have studied uniformly rotating configurations with a poltropic EOS and two tabulated 
"realistic" EOS from the literature. Physical characteristics of the maximum mass configurations 
for both nonrotating and rotating stars are listed in Table 3, and compared with the results of 
two other groups. While there is excellent agreement across the board, the agreement is particu- 
larly good in the polytropic case. Slightly larger differences in the case of the tabulated EOS are 
well-attested in the literature, resulting from different methods of interpolation, matching between 
different density regimes, etc. We tried two methods of determining the "heat function" h from the 
tabulated EOS: 1) direct integration of equation (37), and 2) use of an analytic formula derived 
from equation (37) with the first law of thermodynamics. Bocquet et al. (1995) used the analytic 
formula, and when employing this method we obtained closer agreement with the results of those 
authors. However, we achieved smaller values of |1 — A| (indicating better solutions) when con- 
structing h by direct integration. This was apparently the method employed by Cook, Shapiro, & 
Teukolsky (1994b), as our results with this method agree more closely with theirs. Our calcula- 
tions underlying the entries in Table 3 — and, indeed, our calculations with tabulated EOS reported 
throughout this work — employed the construction of h by direct integration. 

In our experience the calculations of rotating configurations with tabulated EOS had a tendency 
towards numerical instability for large angular velocities, manifested as growing oscillations in 
various quantities (as opposed to the monotonic runaway occurring past the Keplerian limit). We 
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found that this instability could be controlled by updating the metric variable N'f' only once every 
n iterations, where we took n as high as 15 when approaching the Keplerian limit. This procedure 
was not necessary with the polytropic EOS we used. 

We also tested our code's ability to reliably compute static neutron star configurations with 
large magnetic fields. As reported in the main text, our results for the maximum mass of neutron 
stars with constant current function /o are larger than those reported by Bocquet et al. (1995). 
Hence comparison of these maximum mass configurations does not really constitute a verification of 
our code. However, Bocquet et al. (1995) also presented results for the maximum mass at a certain 
fixed low values of magnetic moment, and comparison of these results provides a benchmark against 
which our code can be checked. In addition, we can compute configurations close to those reported 
by Bocquet et al. (1995) as the maximum mass, and make a comparison — even though these are 
not our maximum mass configurations. 

The results of these calculations are presented in Table 4. For each EOS, three sets of calcula- 
tions arc presented. The first set contains results of stars without magnetic fields — the same data 
appearing in Table 3 — as a baseline. The second set shows results for the maximum mass at a given 
fixed (relatively low) value of magnetic moment. The third set for each EOS has two entries of our 
calculations. These two configurations represent the boundaries of the range of configurations, at 
fixed magnetic moment, whose baryon mass rounds to the value reported by Bocquet et al. (1995) 
as the maximum baryon mass among all configurations with constant current function. The results 
are satisfactorily concordant. For a visual comparison, our Figure 1 can be checked against Figures 
5 and 6 of Bocquet et al. (1995). 

We have cited three classes of tests which validate our code. First, the quantity |1 — A| is close 
to zero as required of valid solutions. Second, our results for quantities characterizing the maximum 
mass configuration among uniformly rotating stars agree well with those of two previous groups. 
Third, quantities characterizing certain configurations with magnetic fields reported by Bocquet et 
al. (1995) show good agreement. 
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Table 1. Properties of soft 7-ray repeaters (SGRs) and anomalous X-ray pulsars (AXPs) from 
the recent literature*. Question marks indicate uncertain or unconfirmed values. 



SGR 


P 


p a 


B b 


D " 




Supernova 


Comments 




(s) 


(10-" s/s) 


(10" G) 


(kpc) 


(1034 erg/s) 


Remnant (SNR) 




1806-20 


7.48 


2 


8 


14 


20 


GlO.0-0.3 




1900+14 


5.16 


6.1 


8 


5 


3 


G42.8+0.6? 


27-Aug-98 giant flare; 
















radio point source 


0525-66 


8 




4.5 


55 


~ 300 


N49 in LMC 


05-Mar-79 giant flare 


1627-41 


6.4? 






11 


10 


G337.0-0.1 




1801-23 








10? 






very recent; only two bursts 


AXP 


4U 0412+61 


8.69 


0.23 


2.8 


4 


112 


none 


Lbb ~ 40% Lx- P constant. 


IE 1048.1-5937 


6.45 


1.67/3.29/1.67 


6.5 


10 


8 


none 


^BB ~ 55% Lx- Three epochs 
















of different, but constant, P. 


IRXS J170849-4009 


11.00 


1.9 


9.2 


10 


140 


none 


-f'BB 55% Lx- Regular 
















spindown, except for a glitch. 


IE 1841-045 


11.77 


4.1 


14 


7 


35 


Kes 73 




AX J1845.0-0258 


6.97 






15 


30 


G29.6+0.1 


No P to date but young SNR, 
















hence large P. 


IE 2259+586 


6.97 


0.06 


1.3 


4 


3.3 


CTB 109 


Lbb ~ 40% Lx- Bumpy 
















spin-down. 


Magnctar candidates 


RX J0720.4-3125 


8.39 






0.1 


0.0026 


none 


Blackbody (BB) spectrum, 
















proposed old magnetar 


RX J0420.0-5022 


22.7? 






4 


0.4 


none 


Needs confirmation. 
















Candidate because of large P. 



*The entries in this table axe extracted from Cline, et al. (2000); Corbel, Chapuis, Dame, & Durouchoux (1999); Gaensler, Gotthelf, & 
Vasisht (1999); Gotthelf & Vasisht (1998); Gotthelf, Vasisht, & Dotani (1999); Haberl, Motch, Buckley, Zickraf, & Pietsch (1997); Habcrl, 
Pietsch, & Motch (2000); Hurley (2000); Hurley et al. (2000); Heyl & Hcrnquist (1998); Kaspi, Lackey, & Chakrabarty (2000); Kouvcliotou 
et al. (1998, 1999); Oosterbroek, Parmar, Mercghetti, & Israel (1998); Parmar et al. (1998); Paul, Kawasaki, Dotani, & Nagase (2000); 
Rothschild, Kulkarni, & Lingenfelter (1994); Sugizaki, et al. (1997); Torii et al. (1998); Vasisht & Gotthelf (1997); Vasisht, Kulkarni, Anderson, 
& Hamilton (1997); White et al. (1996); Woods et al. (1999a,b) 

^'In most cases there is not enough sampling of P to allow strong claims of constancy in time, and in many cases there is evidence of 
significant variations. 

is obtained by standard magnetodipolar radiation braJcing, only good to within an order of magnitude. 

■^Distances are poorly determined (with the exception of SGR 0525-66 in the LMC). 

''The X-ray luminosity Lx scales as D"^ and is estimated from the observed unabsorbed flux using the quoted distance. 
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Table 2. Maximum mass models at various fixed values of magnetic moment. 



EOS 


M 


^max 


Gmax 


Be 


-^polc 


M 


Mb 


R 


|1-A| 




(10^^ Gaussian) 




(1.66 X 1014 g cm-3) 


(IOI6 G) 


(IOI6 G) 


(M0) 


(Mo) 


(km) 




rOiZ 


U.UU 


n /loo 


lU.O 


U.UU 


U.UU 


nn 
z.UU 


010 
z.iy 


1 Q Q 

io.cS 


IT? A 




2.00 


0.459 


9.41 


150 


36.9 


2.12 


2.28 


13.9 


6E-5 




3.00 


0.511 


11.1 


253 


122 


2.37 


2.38 


12.8 


2E-4 




4.00 


0.386 


7.28 


180 


92.9 


2.56 


2.66 


15.2 


3E-4 




4-99 


0.291 


4.94 


129 


70.6 


2.64 


2. 78 


17.6 


3E-4 




5.00 


0.290 


4.92 


129 


7U.4 


2.64 


2.78 


1 T ft 

17.6 


3E-4 




6.18 


0.143 


2. 07 


59.6 


34.4 


2. 31 


2.43 


22.9 


3E-4 


rSJl 


U.UU 


U.D<5 1 


1 Q C 

icS.O 


n nn 
U.UU 


n nn 
U.UU 


i.OO 




n 00 


A 




1.50 


0.568 


14.6 


236 


84.1 


1.96 


2.17 


10.3 


4E-5 




2.00 


0.706 


19.2 


419 


249 


2.12 


2.05 


9.53 


AT? A 




2.50 


0.586 


15.1 


335 


216 


2.28 


2.29 


10.7 


5E-4 




3.00 


0.476 


11.9 


269 


181 


2.38 


2.47 


11.9 


6E-4 




3.29 


0.40 8 


10.1 


232 


157 


2.40 


2.54 


12.6 


lE-3 




3.70 


0. 265 


6.68 


156 


104 


2.25 


2.41 


14.4 


2E-4 


FanclIN 


0.00 


0.728 


24.9 


0.00 


0.00 


1.66 


1.92 


8.36 


5E-5 




1.00 


0.617 


20.3 


292 


87.8 


1.72 


1.92 


8.72 


8E-5 




1.50 


0.724 


24.7 


504 


303 


1.86 


1.80 


8.18 


3E-4 




2.00 


0.590 


19.3 


386 


271 


2.05 


2.06 


9.37 


4E-4 




2.49 


0.422 


13.7 


281 


202 


2.13 


2.23 


10.8 


5E-4 




2.50 


0.410 


13.4 


276 


197 


2.12 


2.23 


10.8 


5E-4 




2.82 


0.231 


8. 31 


166 


119 


1.90 


2.03 


12.6 


2E-4 


Akmal 


n An 
U.UU 


u.yiu 


10. / 


n nn 
U.UU 


n nn 
U.UU 


on 
Z.ZKJ 


ft'? 


1 n n 
iU.U 


IT? A 




1.20 


0.792 


14.3 


228 


59.3 


2.22 


2.62 


10.3 


7E-5 




2.00 


0.648 


11.8 


338 


134 


2.31 


2.58 


10.6 


3E-4 




2.80 


0.674 


12.3 


370 


244 


2.53 


2.59 


11.0 


2E-4 




3.60 


0.550 


10.3 


287 


221 


2.73 


2.84 


12.4 






3.71 


0. 514 


9.75 


271 


209 


2.73 


2.88 


12.6 


3E-4 




4.10 


0. 307 


6. 81 


179 


132 


2.52 


2.73 


14.3 


lE-3 


PCL 


0.00 


0.561 


17.3 


0.00 


0.00 


1.72 


1.95 


10.4 


4E-4 




1.70 


0.665 


22.7 


418 


218 


1.91 


1.86 


9.14 


4E-4 




2.20 


0.547 


16.6 


327 


186 


2.09 


2.12 


10.4 


6E-4 




2.70 


0.445 


12.5 


259 


155 


2.21 


2.31 


11.7 


8E-4 




3.20 


0.350 


9.22 


201 


124 


2.26 


2.41 


13.1 


5E-4 




3.31 


0.335 


8.75 


192 


120 


2.28 


2.43 


13.4 


4E-4 




3.86 


0.202 


5.12 


116 


74.7 


2.09 


2.24 


15.8 


IE- 4 


PCLhyp 


0.00 


0.504 


17.6 


0.00 


0.00 


1.59 


1.78 


10.3 


4E-4 




1.50 


0.620 


24.3 


416 


202 


1.76 


1.73 


8.82 


5E-4 




2.00 


0.488 


16.7 


308 


161 


1.93 


1.98 


10.3 


5E-4 




2.50 


0.387 


12.1 


235 


130 


2.04 


2.15 


11.8 


3E-4 




3.00 


0.310 


8.97 


182 


106 


2.11 


2.25 


13.2 


8E-6 




3.39 


0.263 


7.29 


152 


93.1 


2.13 


2.28 


14.2 


2E-4 




3.81 


0.186 


4.83 


107 


68.7 


2.02 


2.16 


16.1 


5E-4 


Note. — 


The quantity hmax 


is the 


maximum value of the log-cnthalpy per baryon 


defined in 


equation 


(37), 


^max is the 



maximum energy density, Ba and Bp^ie are the magnetic field strengths at the star's center and pole, respectively, M is the 
gravitational mass, Mb is the star's baryon mass, R is the equatorial radius, and |1 — A| is a function defined in equation 
(A13) describing convergence. The results are grouped by the equation of state (EOS), each of which is described in the text. 
Entries in plain type correspond to the sequences of constant magnetic moment M shown in Figures 7 and 8. For each EOS 
there are two italicized rows which correspond, respectively, to the configurations of maximum mass and maximum magnetic 
moment among all configurations with constant current function. 
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Table 3. Maximum mass models, nonrotating and rotating. 



EOS 


Authors 


he 


ec 


n 


M 


Mb 


R 


J/M2 


|1-A| 








(1.66 X 10" g cm-3) 


(10* s-1) 


(Mq) 


(Mq) 


(km) 


(Gjv/c) 




Pol2 


GST 




10.6 


0.00 


1.99 


2.19 


13.7 


0.00 






BBGN 


0.491 


10.4 


0.00 


2.00 


2.19 


13.8 


0.00 


IE- 14 




CPL 


0.492 


10.5 


0.00 


2.00 


2.19 


13.8 


0.00 


lE-4 




CST 




8.63 


0.629 


2.29 


2.52 


19.7 


0.572 






BBGN 


0.432 


8.58 


0.629 


2.30 


2.52 


19.6 


0.570 


6E-6 




GPL 


0.432 


8.58 


0.629 


2.. 30 


2.53 


19.6 


0.570 


lE-4 


BJI 


GST 




18.5 


0.00 


1.86 


2.16 


9.93 


0.00 






BBGN 


0.699 


18.6 


0.00 


1.86 


2.13 


9.91 


0.00 


2E-6 




GPL 


0.687 


18.5 


0.00 


1.86 


2.14 


9.92 


0.00 


5E-4 




GST 




16.0 


1.06 


2.17 


2.49 


13.4 


0.629 






BBGN 


0.628 


16.3 


1.07 


2.15 


2.46 


13.4 


0.626 


9E-5 




GPL 


0.610 


15.9 


1.07 


2.16 


2.47 


13.4 


0.632 


4E-5 


PandN 


GST 




24.9 


0.00 


1.66 


1.92 


8.37 


0.00 






BBGN 


0.733 


24.4 


0.00 


1.66 


1.93 


8.55 


0.00 


2E-4 




GPL 


0.728 


24.9 


0.00 


1.66 


1.92 


8.36 


0.00 


5E-5 




CST 




21.4 


1.32 


1.95 


2.24 


11.2 


0.665 






BBGN 


0.668 


21.7 


1.29 


1.93 


2.23 


11.4 


0.641 


7E-5 




GPL 


0.647 


21.5 


1.33 


1.96 


2.25 


11.1 


0.666 


lE-5 



Note. — J is the total angular momentum. See the caption to Table 2 for an explanation of the other quantities 
tabulated. The subscript c refers to central values. For each EOS, results from this work (GPL) are compared 
to those of Gook, Shapiro, & Teukolsky (1994a) and Gook, Shapiro, & Teukolsky (1994b) for the polytropic and 
tabulated EOSs, respectively (together labelled GST), and to Bocquet et al. (1995) (BBGN). 
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Table 4. Comparison of static models with various magnetic moments. 



EOS 


Authors 


M 


he 


ec 


Be 


Bpole 


M 


Mb 


R 


|1-A| 






(10^5 Gaussian) 




(1.66 X 10" g cm-3) 


(10i« G) 


(10i« G) 


(Mo) 


(Mo) 


(km) 




Pol2 


CST 


0.00 




10.6 


0.00 


0.00 


1.99 


2.19 


13.7 






BBGN 


0.00 


0.491 


10.4 


0.00 


0.00 


2.00 


2.19 


13.8 


lE-14 




CPL 


0.00 


0.492 


10.5 


0.00 


0.00 


2.00 


2.19 


13.8 


lE-4 




BBGN 


0.800 


0.483 


10.2 


59.8 


9.8 


2.01 


2.21 


13.8 


lE-6 




CPL 


0.800 


0.484 


10.2 


60.4 


10.0 


2.02 


2.21 


13.8 


2E-4 




BBGN 


4.49 


0.225 


3.55 


142 


72.3 


2.57 


2.71 


16.8 


lE-3 




CPL 


4.49 


0.224 


3.53 


140 


70.7 


2.56 


2.71 


16.8 


2E-4 




CPL 


4.49 


0.224 


3.53 


142 


72.6 


2.57 


2.71 


16.8 


3E-4 


BJI 


CST 


0.00 




18.5 


0.00 


0.00 


1.86 


2.16 


9.93 






BBGN 


0.00 


0.699 


18.6 


0.00 


0.00 


1.86 


2.13 


9.91 


2E-6 




CPL 


0.00 


0.687 


18.5 


0.00 


0.00 


1.86 


2.14 


9.92 


5E-4 




BBGN 


0.30 


0.692 


18.4 


61.5 


11.0 


1.86 


2.13 


9.94 


2E-6 




CPL 


0.30 


0.680 


18.2 


60.9 


11.1 


1.86 


2.14 


9.95 


2E-4 




BBGN 


2.63 


0.300 


7.47 


233 


121 


2.18 


2.34 


12.0 


lE-4 




CPL 


2.63 


0.287 


7.16 


215 


110 


2.15 


2.34 


12.3 


2E-4 




CPL 


2.63 


0.292 


7.27 


222 


115 


2.17 


2.34 


12.2 


2E-4 


PandN 


CST 


0.00 




24.9 


0.00 


0.00 


1.66 


1.92 


8.37 






BBGN 


0.00 


0.733 


24.4 


0.00 


0.00 


1.66 


1.93 


8.55 


lE-4 




CPL 


0.00 


0.728 


24.9 


0.00 


0.00 


1.66 


1.92 


8.36 


5E-5 




BBGN 


0.20 


0.727 


24.1 


64.7 


11.9 


1.66 


1.93 


8.57 


lE-4 




CPL 


0.20 


0.722 


24.6 


66.6 


12.9 


1.66 


1.92 


8.38 


4E-5 




BBGN 


1.86 


0.350 


11.2 


303 


154 


1.91 


2.06 


10.0 


3E-4 




CPL 


1.86 


0.328 


11.0 


292 


153 


1.90 


2.06 


10.0 


3E-4 




CPL 


1.86 


0.336 


11.2 


361 


209 


1.96 


2.06 


9.47 


3E-4 



Note. — 



See the notes to Tables 2 and 3 for explanations of symbols and abbreviations. 
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Fig. 1. — Contour plots of the energy density and the electromagnetic vector potential component 
Aff). Here x = rsinO and y = rcosO, where r and 6 are coordinates appearing in equation (2); note 
that in these coordinates not even distances in the equatorial plane constitute proper distances, 
eo = 1.66 X 10^^ g cm~^, and i?* is the equatorial radius. While contours of constant show 
the structure of the magnetic field, their spacing as pictured here does not accurately indicate 
magnetic field strength; the maximum magnetic field is actually at the center of the star. To allow 
a direct comparison, with Figures 5 and 6 of Bocquct ct al. (1995), we assumed their value of the 
polytropic constant. Table 4 contains physical quantities rescaled to reflect a more realistic value 
of this constant). 
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Fig. 2. — Various quantities associated with the Bcrnoulh equation (42) in the equatorial plane of 
a uniformly rotating star, as a function of the compactified radial coordinate s = r/{R^, + r), where 
R* is the equatorial radius. The upper panel shows a configuration with modest rotation, and the 
lower panel shows a configuration at the Keplerian limit. See the text for discussion. 
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Electromagnetic Current Forces 




Fig. 3. — The electromagnetic current and various "forces" in the equatorial plane for the config- 
uration pictured in Figure 1. The "forces" are derivatives of the terms in the Bernoulli equation 
(43); here M is the magnetic potential. The right boundaries correspond to the equatorial radius 
i?*. Note that the pressure force {—dh/dr) points inward at lower radii; this is a consequence of the 
maximum density being pushed off center (see the left panel of Figure 1) due to the strong outward 
Lorentz force {—dM/dr). The Lorentz force also reverses sign, due to a reversal of magnetic field 
direction in the equatorial plane (see the right panel of Figure 1). At i?*, the Lorentz force works 
together with gravity to help confine the star, in contrast with the centrifugal force in the case of 
rapid rotation. 
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Fig. 4. — Similar to Figure 2, but for the Bernoulli equation (43) of a nonrotating star in a poloidal 
magnetic field. M is the magnetic potential. The upper panel shows a configuration in which the 
magnetic field is strong enough to modestly deform the star, while the lower panel displays the 
configuration with the largest current function for which convergence was achieved (for the same 
maximum density as the configuration in the upper panel). When compared with Figure 2, the 
absence of mass shedding is evident. 
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Fig. 5. — Late stages of iteration of a non-converging configuration (as a solution of the Einstein 
and Maxwell equations). Left panels: Density contour plots show a transition to a toroidal topology 
as the iterations proceed. Center panels: "Potential" plots similar to Figure 2; note the increasingly 
narrow well into which the matter is compressed. Right panels: "Force" plots similar to Figure 4; 
the tendency of the gravitational force to join the magnetic force in pushing matter away from the 
origin becomes more and more pronounced. 
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Fig. 6. — The gravitational force in the equatorial plane. Upper panel: The gravitational force for 
increasing values of /o, at a fixed value of maximum stellar density. The indicated values of /o 
are given in units of 10^^(i?* eoc^)~^ A m~^, where is the coordinate radius of the equatorial 

surface. The insets show the force near the origin, computed both "analytically" and numerically, 
as described in the text. Lower panel: The gravitational force for the maximum value of /o for 
which convergence was achieved (for a particular value of maximum density), for configurations 
computed with three different resolutions in r and 9. 
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Fig. 7. — Mass-equatorial radius plots showing converged solutions attainable with a constant 
current function, for EOSs used by Bocquet et al. (1995). The lower heavy line represents spherical, 
non-magnetized, configurations, and the upper heavy line represents the boundary beyond which 
solutions appear not to exist (see §3). Lighter solid lines are sequences of constant baryon mass (in 
Mq), while dotted lines are sequences of constant M. (in units of Al* = 10^^ Gaussian). Lighter 
shaded regions indicate configurations in which the maximum density is not at the center; the 
darker shaded regions indicate gravitationally unbound configurations. "X"s denote the maximum 
mass configuration attainable by uniform rotation, and crosses indicate the maximum masses for 
non-rotating, magnetized configurations reported by Bocquet et al. (1995). 
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Fig. 8. — Similar to Figure 7, but for three more recent EOSs discussed in the text. 



